using DataFrames
using QuadGK
using Distributions
using CSV
using LinearAlgebra
using JLD
using DelimitedFiles
using SparseArrays

clearconsole()

#-------------------------------------------------------------
# IRF Configuration
#-------------------------------------------------------------

H = 40 # maximum horizon
n_every = 10 # use every n_every'th draw from the posterior
sh_size = 3   # shock size, in multiples of standard deviations
sh_id = 3 # sh_id = 1: TFP shock

#-------------------------------------------------------------
# include Functions
#-------------------------------------------------------------
#cd("$(pwd())/Dropbox/Heterogeneity/Software/KS_Simulation/")
readDir = "$(pwd())/Functions/"
# include(readDir *"vech.jl");
include(readDir *"logSpline_Procedures.jl");
include(readDir *"VAR_MoreLags_Procedures.jl");
include(readDir *"Loaddata.jl");

#-------------------------------------------------------------
# choose specification files
#-------------------------------------------------------------
nVARSpec    = "1"
nVARMCMCSpec= "1"
specDir   = "$(pwd())/SpecFiles/"
include(specDir * "/AggVARspec" * nVARSpec * ".jl")
include(specDir * "/AggVARMCMCspec" * nVARMCMCSpec * ".jl")

#-------------------------------------------------------------
# load aggregate data
#-------------------------------------------------------------
dataDir  = "$(pwd())/data/"
agg_data, period_agg, ~ = loadaggdata(SampleStart,SampleEnd)
n_agg    = size(agg_data)[2]
data_adj = agg_data[Tdrop-nlags+1:end,:]

n = n_agg
p = nlags


#-------------------------------------------------------------
# Load posterior draws
#-------------------------------------------------------------
sName    = "AggVAR" * nVARSpec * "_MCMC" * nVARMCMCSpec
loadDir  = "$(pwd())/results/" * sName *"/";

Bpdraw     = load(loadDir * sName * "_PostDraws.jld", "Bpdraw")
Sigtrpdraw = load(loadDir * sName * "_PostDraws.jld", "Sigtrpdraw")
Bpmean     = load(loadDir * sName * "_PostMeans.jld", "Bpmean")
Sigtrpmean = load(loadDir * sName * "_PostMeans.jld", "Sigtrpmean")

#-------------------------------------------------------------
# Generate IRFs at posterior mean of (Phi,Sigmatr)
#-------------------------------------------------------------

println("")
println("Generating IRFs at posterior mean... ")
println("")

YY_IRF   = zeros(H+1,n)
L        = Sigtrpmean
Btilde   = reshape(Bpmean, n*p,n)
response = IRFredu(Btilde,L,H,sh_size,sh_id)
YY_IRF[2:end,:] = response'

savedir = "$(pwd())/results/" * sName *"/";
try mkdir(savedir) catch; end
CSV.write(savedir * sName * "_IRF_YY_AggSh"*string(sh_id)*"_pmean.csv", DataFrame(YY_IRF,:auto))

#-------------------------------------------------------------
# Generate IRFs for a subset of posterior draws
#-------------------------------------------------------------

println("")
println("Generating IRFs for posterior draws... ")
println("")

nsim = size(Bpdraw)[1]
n_subseq                 = floor(Int,nsim/n_every)
YY_IRF_uncertainty       = zeros(H+1,n,n_subseq)

# hh=1 is steady state
# hh=2 is period of impact

for pp = 1:n_subseq

    time_init_loop = time_ns();
    println(" Draw number:  $pp")
    println(" Remaining draws:  $(n_subseq-pp)")
    println(" Posterior draw number: $(pp*n_every)")

    L        = Sigtrpdraw[pp*n_every,:,:]
    Btilde   = Bpdraw[pp*n_every,:,:]
    response = IRFredu(Btilde,L,H,sh_size,sh_id)
    YY_IRF_uncertainty[2:end,:,pp] = response'

    CSV.write(savedir * sName * "_IRF_YY_AggSh" * string(sh_id) * "_" * string(pp) * ".csv", DataFrame(YY_IRF_uncertainty[:,:,pp],:auto))

    time_loop=signed(time_ns()-time_init_loop)/1000000000
    println(" Elapsed time = $(time_loop)")
    println("")

end
